O.4 Frontiers OXT - Review

Updates for Review

1 Read and Merge

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
wd <- "C:/Users/Mihai/Desktop/O.4 prealabil pt Frontiers/O.4 Date Psiho&OXT"
setwd(wd)

Data_OXT <- rio::import(file.path(wd, "\\OXT", "O.4 OXT Data.xlsx"))

Data_psih <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "26martie2019")
Data_psih <- Data_psih[, c(1:6, 8:13)]

Data_Neo <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "Neo scored")
Data_BDI <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "BDI scored")
Data_STAI <- rio::import(file.path(wd, "\\Psiho", "DATE O4Cl 16.11.2019 procesat.xlsx"), which = "STAI scored")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Clean & Merge
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Data_STAI <- 
  Data_STAI %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 2) %>%                                                     # filter out rows (no more than 1 NA on row)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_BDI <- 
  Data_BDI %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 1) %>%                                                     # filter out rows (not 1 NA)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_Neo <- 
  Data_Neo %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 6) %>%                                                     # filter out rows (no more than 5 NA on row)
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::select(-Experim)

Data_psih <- 
  Data_psih %>% 
  dplyr::filter(!is.na(ID)) %>%                                                                # no ID, no data
  dplyr::filter(rowSums(is.na(.)) < 8) %>%                                                     # filter out rows (no more than 7 NA on row)
  tidyr::separate(Conditia,  c("Nr_zi", "Conditia"), "\\s+") %>%                               # split on white space
  mutate(ID = stringr::str_replace(ID, "\\s", "|")) %>%                                        # replace first whitespace with | and use it to separate
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\|") %>%                              # separate on only first whitespace that was replaced
  dplyr::rename("Nume" = Nume_Prenume) %>% 
  dplyr::select(-Experim)

Data_OXT <-   
  Data_OXT %>% 
  tidyr::separate(ID,  c("ID", "Ziua", "Nr_zi", "Proba"), "\\s+") %>%                          # split on white space
  select(ID, Nr_zi, Proba, OXT) %>%
  tidyr::separate(ID, into = c("ID", "Experim"), sep = "\\/") %>%                              # separate on /
  dplyr::select(-Experim) %>% 
  dplyr::rename("PrePost" = Proba) %>% 
  mutate(PrePost = forcats::fct_recode(PrePost, "Pre" = "A", "Post" = "B")) %>%                # Proba A = Pre, B = Post on same day
  spread(key = PrePost, value = OXT) %>%
  dplyr::rename_at(vars(Pre, Post), ~ c("OXT_Pre", "OXT_Post"))


Data_merged <- dplyr::left_join(Data_psih, Data_OXT, by = c("ID", "Nr_zi"))                 # ID 40, 59 are incomplete in OXT -- dplyr::full_join


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Make Wide Dataframe of merged
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Check if ids have > 1 row of data (empty ID have only 1 row; thei have only 1 Condition) 
# Careful! This function modfies the datasets in the global envinronment
delete_empty_id <- function(df){
  list_empty_id <- 
    df %>%
    dplyr::group_by(ID) %>%
    dplyr::summarise(row_count = n()) %>%
    dplyr::rename("empty_id" = ID) %>%
    mutate(delete_id = if_else(row_count < 2, TRUE, FALSE)) %>%
    filter(delete_id == TRUE)
  
  df_modif <- 
    df %>%
    filter(!ID %in% list_empty_id$empty_id)
  
  if(!identical(df, df_modif)){
    df <- deparse(substitute(df))
    cat("Deleting from ", print(as.name(df))); print(list_empty_id)                    # print out which ids are deleted from which dataset
    assign(df, df_modif, envir = globalenv())                                          # assign modified df to original dataset from Global
  }else cat("No empty datasets. Nothing to delete")
}
delete_empty_id(Data_merged) 

Data_merged Deleting from Data_merged

# Dataframe for regular analyses
Data_merged_wide <-
  Data_merged %>%
  select(-c(Varsta, Gen)) %>%            # exclude Gen, Varsta: merge after because it interferes with spread()
  gather("variable", "value", c(Data, Nr_zi, Vas_Stres_Pre:OXT_Post), -c(1:2)) %>%       # Conditia needs to be outside
  unite("united_var", c(variable, Conditia), sep = "_") %>%
  spread(united_var, value) %>%
  mutate_at(vars(5:22), list(~as.numeric(as.character(.))))

# sum(is.na(Data_merged[, 6:15])) == sum(is.na(Data_merged_wide[, 5:22]))    # check if nr of NA is the same after melt


# Full Wide Dataframe  (HERE We APPLY EXCLUSION CRITERIA) 
Data_merged_wide_all <- 
  Data_merged[, c("ID", "Nume", "Varsta", "Gen")] %>%                        # need to add these back, but just the main unduplicated row
  dplyr::distinct(ID, .keep_all = TRUE) %>%
    dplyr::left_join(., Data_merged_wide, by = "ID") %>%
    dplyr::rename("Nume" = Nume.x) %>%
    dplyr::select(-Nume.y) %>% 
    dplyr::left_join(., Data_BDI, by = "ID") %>%
    dplyr::left_join(., Data_STAI, by = "ID") %>%
    dplyr::left_join(., Data_Neo, by = "ID") %>%
    mutate_at(vars(23:60), list(~as.numeric(as.character(.)))) %>%
    mutate(Varsta = as.numeric(as.character(Varsta))) %>% 
    filter(!ID %in% c(19, 43)) %>%           # 55 subs remain after this
    filter(Varsta < 40) %>%                  # 53 subs remain after this
    filter(ScorBDI < 30) %>%                 # 50 subs remain after this
    dplyr::arrange(ID)


# Full Long Dataframe
vars_OglEcran <- grep("OGL|ECRAN", colnames(Data_merged_wide_all))

Data_merged_long_all <-
  Data_merged_wide_all %>%
  gather(variable, value, vars_OglEcran, -c(ID)) %>%                                      # gather Data and Nr_zi as these are OglEcran level
  tidyr::separate(variable,  c("variable", "OglEcran"), "_(?=[^_]+$)") %>%                # split only on last "_"
  spread(variable, value) %>% 
  rename_at(vars(vars_OglEcran), list(~stringr::str_remove_all(., c("_OGL|_ECRAN")))) %>%          
  arrange(ID)

vars_PrePost <- grep("Pre|Post", colnames(Data_merged_long_all))

Data_merged_long_all <-
  Data_merged_long_all %>%
  pivot_longer(vars_PrePost, names_to = c(".value","PrePost"), names_pattern = "(.*)_(.*)") %>%     # magic     
  mutate_at(vars(48:50), list(~as.numeric(as.character(.)))) %>%
  arrange(ID) %>%
  rstatix::convert_as_factor(ID, PrePost, OglEcran)


# Females and Males Dataframes
Data_fem <- 
  Data_merged_long_all %>%
  filter(Gen == "f")
Data_masc <- 
  Data_merged_long_all %>%
  filter(Gen == "m")

# Difference Scores Dataframe
Data_difscores <-
  Data_merged_wide_all %>%
  mutate(Diff_OXT_OGL = OXT_Post_OGL - OXT_Pre_OGL,
         Diff_OXT_ECRAN = OXT_Post_ECRAN - OXT_Pre_ECRAN,
         Diff_Vas_Stres_OGL = Vas_Stres_Post_OGL - Vas_Stres_Pre_OGL,
         Diff_Vas_Stres_ECRAN = Vas_Stres_Post_ECRAN - Vas_Stres_Pre_ECRAN, 
         Diff_Vas_Bine_OGL = Vas_Bine_Post_OGL - Vas_Bine_Pre_OGL,
         Diff_Vas_Bine_ECRAN = Vas_Bine_Post_ECRAN - Vas_Bine_Pre_ECRAN,
         Diff_IOS_OGL = IOS_Post_OGL - IOS_Pre_OGL,
         Diff_IOS_ECRAN = IOS_Post_ECRAN - IOS_Pre_ECRAN
  ) 

1.1 Dataset - Long Format

Data_merged_long_all %>%
  dplyr::select(-Nume) %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

1.2 Dataset - Wide Format

Data_merged_wide_all %>%
  dplyr::select(-Nume) %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

2 Define functions

# Define Function for Two-way rmANOVA
# library(tidyverse)
# library(ggpubr)
# library(rstatix)
# library(rlang)

tw_rmANOVA_func <- 
  function(data, id_var, cond_var, time_var, value_var, 
           assum_check = TRUE, posthoc_sig_interac = FALSE, posthoc_ns_interac = FALSE,
           p_adjust_method = "bonferroni"){
  
  # input dataframe needs to have columns names diffrent from "variable" and "value" because it collides with rstatix::shapiro_test
    
  id_var_enq <- rlang::enquo(id_var)  
  cond_var_enq <- rlang::enquo(cond_var)
  cond_var_name <- rlang::as_name(cond_var_enq)
  time_var_enq <- rlang::enquo(time_var)
  time_var_name <- rlang::as_name(time_var_enq)
  value_var_enq <- rlang::enquo(value_var)
  value_var_name <- rlang::as_name(value_var_enq)
    
  # Assumptions
  if(assum_check){
    cat("\n Outliers \n")
    data %>%
      dplyr::group_by(!!cond_var_enq, !!time_var_enq) %>%
      rstatix::identify_outliers(!!value_var_enq) %>%                                  # outliers (needs to be 0)
      print()
    
    cat("\n Normality assumption (p>.05) \n")
    data %>%
      dplyr::group_by(!!cond_var_enq, !!time_var_enq) %>%
      rstatix::shapiro_test(!!value_var_enq) %>%                                        # normality assumption (p>.05)
      print()
    
    qq_plot <- 
      ggpubr::ggqqplot(data = data, value_var_name, ggtheme = theme_bw(), title = "QQ Plot") +
        ggplot2::facet_grid(vars(!!time_var_enq), vars(!!cond_var_enq), labeller = "label_both")    # QQ plot
  }
  
  # Two-way rmANOVA - check for interaction (ex. F(2, 22) = 30.4, p < 0.0001)
  cat("\n Two-way rmANOVA \n")
  res_aov <- anova_test(                         # automatically does sphericity Mauchly’s test
    data = data, dv = !!value_var_enq, wid = !!id_var_enq,                             
    within = c(!!cond_var_enq, !!time_var_enq)
  )
  get_anova_table(res_aov) %>%  # ges: Greenhouse-Geisser sphericity correction is automatically applied to factors violating the sphericity assumption  
  print()
  
  #- Procedure for a significant two-way interaction -
  if(posthoc_sig_interac){
    cat("\n Effect of treatment at each time point \n")
    one_way <- 
      data %>%
      group_by(!!time_var_enq) %>%
      anova_test(dv = !!value_var_enq, wid = !!id_var_enq, within = !!cond_var_enq) %>%
      get_anova_table() %>%
      adjust_pvalue(method = "bonferroni")
    one_way %>% print()
    
    cat("\n Pairwise comparisons between treatment groups \n")
    pwc <-
      data %>%
      group_by(!!time_var_enq) %>%
      pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", cond_var_name)),
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc %>% print()

    cat("\n Effect of time at each level of treatment - One-way ANOVA \n")
    one_way2 <-
      data %>%
      group_by(!!cond_var_enq) %>%
      anova_test(dv = !!value_var_enq, wid = !!id_var_enq, within = !!time_var_enq) %>%
      get_anova_table() %>%
      adjust_pvalue(method = p_adjust_method)
    one_way2  %>% print()
    cat("\n Pairwise comparisons between time points \n")
    pwc2 <-
      data %>%
      group_by(!!cond_var_enq) %>%
      pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", time_var_name)),     # paste formula, not quosure
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc2  %>% print()
  }
  
  #- Procedure for non-significant two-way interaction- 
  # If the interaction is not significant, you need to interpret the main effects for each of the two variables: treatment and time.
  if(posthoc_ns_interac){
    cat("\n Comparisons for treatment variable \n")
    pwc_cond <-
      data %>%
      pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", cond_var_name)),     # paste formula, not quosure             
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc_cond %>% print()
    cat("\n Comparisons for time variable \n")
    pwc_time <-
      data %>% 
      pairwise_t_test(
        as.formula(paste0(value_var_name, " ~ ", time_var_name)),     # paste formula, not quosure
        paired = TRUE,
        p.adjust.method = p_adjust_method
      )
    pwc_time %>% print()
  }
  
  # Visualization
  bx_plot <- 
    ggboxplot(data, x = time_var_name, y = value_var_name,
    color = cond_var_name, palette = "jco")
  pwc <- 
    pwc %>% 
    add_xy_position(x = time_var_name)
  bx_plot <- 
    bx_plot + 
    stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) +
    labs(
      subtitle = get_test_label(res_aov, detailed = TRUE),
      caption = get_pwc_label(pwc)
    )
  
  if(assum_check){ 
    list(qq_plot, bx_plot)
  }else{
    bx_plot
  } 
  
}

# ex. - run on long format
# tw_rmANOVA_func(data = selfesteem2, id_var = id, cond_var = treatment, time_var = time, value_var = scores, 
#                 posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

## Wilcoxon signed rank test on paired samples

# library(tidyverse)
# library(ggpubr)
# library(rstatix)
# library(rlang)

wilcoxon_paired_func <- function(data, id_var, time_var, value_var, 
                                 assum_check = TRUE){
  
  id_var_enq <- rlang::enquo(id_var)  
  id_var_name <- rlang::as_name(id_var_enq)
  time_var_enq <- rlang::enquo(time_var)
  time_var_name <- rlang::as_name(time_var_enq)
  value_var_enq <- rlang::enquo(value_var)
  value_var_name <- rlang::as_name(value_var_enq)
  
  cat("\n Summary statistics \n")
  data %>%
    group_by(!!time_var_enq) %>%
    get_summary_stats(!!value_var_enq, type = c("common")) %>%
    print() 
  
  if(assum_check){  # The test assumes that differences between paired samples should be distributed symmetrically around the median.
  data2 <- 
    data %>% 
    group_by(!!id_var_enq) %>% 
    mutate(differences = !!value_var_enq - lag(!!value_var_enq))

  hist_dif_plot <- gghistogram(data2, x = "differences", y = "..density..", 
                               fill = "steelblue", add_density = TRUE)
  }
  
  cat("\n Wilcoxon signed rank test on paired samples \n")
  stat_test <- 
    data  %>%
    wilcox_test(as.formula(paste0(value_var_name, " ~ ", time_var_name)), 
                paired = TRUE) %>%
    add_significance()
  stat_test %>% print()
  
  cat("\n Effect size \n")
  ef_size <-
    data %>%
    wilcox_effsize(as.formula(paste0(value_var_name, " ~ ", time_var_name)), 
                   paired = TRUE)
  ef_size %>% print()
  
  # Visualization
  bx_plot <- ggpaired(data, x = time_var_name, y = value_var_name, id = id_var_name,
                      ylab = value_var_name, xlab = time_var_name)
  
  stat_test <- 
    stat_test %>% 
    add_xy_position(x = time_var_name)
  
  bx_plot <-
    bx_plot + 
      stat_pvalue_manual(stat_test, tip.length = 0) +
      labs(subtitle = get_test_label(stat_test, detailed = TRUE))
  
  if(assum_check){ 
    list(hist_dif_plot, bx_plot)
  }else{
    bx_plot
  }  
  
}

3 rmANOVA OXT

cat("## All subjects")

## All subjects

Data_merged_long_all %>%
  select(ID, OXT, OglEcran, PrePost) %>%
  drop_na() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = OXT, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd     F     p p<.05      ges
1         OglEcran   1  45 2.449 0.125       0.012000
2          PrePost   1  45 4.436 0.041     * 0.007000
3 OglEcran:PrePost   1  45 0.061 0.806       0.000135

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Females")

## Females

Data_fem %>%
  select(ID, OXT, OglEcran, PrePost) %>%
  drop_na() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = OXT, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd     F     p p<.05   ges
1         OglEcran   1  23 1.845 0.188       0.020
2          PrePost   1  23 4.929 0.037     * 0.017
3 OglEcran:PrePost   1  23 1.400 0.249       0.004

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Males")

## Males

Data_masc %>%
  select(ID, OXT, OglEcran, PrePost) %>%
  drop_na() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = OXT, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd     F     p p<.05      ges
1         OglEcran   1  21 0.581 0.454       0.006000
2          PrePost   1  21 0.169 0.685       0.000471
3 OglEcran:PrePost   1  21 2.078 0.164       0.019000

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

4 rmANOVA Stress

cat("## All subjects")

## All subjects

Data_merged_long_all %>%
  select(ID, Vas_Stres, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%                     # for complete cases
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Stres, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd      F        p p<.05      ges
1         OglEcran   1  48  0.908 0.345000       0.003000
2          PrePost   1  48 12.537 0.000899     * 0.006000
3 OglEcran:PrePost   1  48  0.285 0.596000       0.000182

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Females")

## Females

Data_fem %>%
  select(ID, Vas_Stres, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Stres, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd      F        p p<.05      ges
1         OglEcran   1  25  1.214 0.281000       0.011000
2          PrePost   1  25 14.495 0.000812     * 0.015000
3 OglEcran:PrePost   1  25  0.105 0.748000       0.000127

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Males")

## Males

Data_masc %>%
  select(ID, Vas_Stres, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Stres, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd        F     p p<.05        ges
1         OglEcran   1  22 0.000235 0.988       0.00000128
2          PrePost   1  22 0.732000 0.401       0.00032800
3 OglEcran:PrePost   1  22 0.187000 0.670       0.00028800

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

5 rmANOVA Bine

cat("## All subjects")

## All subjects

Data_merged_long_all %>%
  select(ID, Vas_Bine, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%                     # for complete cases
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Bine, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd     F     p p<.05        ges
1         OglEcran   1  48 0.334 0.566       0.00100000
2          PrePost   1  48 6.113 0.017     * 0.00500000
3 OglEcran:PrePost   1  48 0.013 0.910       0.00000795

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Females")

## Females

Data_fem %>%
  select(ID, Vas_Bine, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Bine, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd        F     p p<.05        ges
1         OglEcran   1  25 0.000406 0.984       0.00000414
2          PrePost   1  25 3.280000 0.082       0.00600000
3 OglEcran:PrePost   1  25 0.074000 0.788       0.00010300

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

cat("## Males")

## Males

Data_masc %>%
  select(ID, Vas_Bine, OglEcran, PrePost) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  tw_rmANOVA_func(id_var = ID, cond_var = OglEcran, time_var = PrePost, value_var = Vas_Bine, 
                  posthoc_sig_interac = TRUE, posthoc_ns_interac = TRUE)

 Outliers 

 Normality assumption (p>.05) 

 Two-way rmANOVA 
ANOVA Table (type III tests)

            Effect DFn DFd     F     p p<.05       ges
1         OglEcran   1  22 1.328 0.262       0.0050000
2          PrePost   1  22 2.712 0.114       0.0040000
3 OglEcran:PrePost   1  22 0.019 0.891       0.0000218

 Effect of treatment at each time point 

 Pairwise comparisons between treatment groups 

 Effect of time at each level of treatment - One-way ANOVA 

 Pairwise comparisons between time points 

 Comparisons for treatment variable 

 Comparisons for time variable 

[[1]]

[[2]]

6 Wilcoxon paired on Order of Conditions - OXT

cat("## Whole Sample - OGL first, ECRAN second")

## Whole Sample - OGL first, ECRAN second

Data_merged_long_all %>%
  select(ID, OXT, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = OXT, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Whole Sample - ECRAN first, OGL second")

## Whole Sample - ECRAN first, OGL second

Data_merged_long_all %>%
  select(ID, OXT, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = OXT, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - OGL first, ECRAN second")

## Females - OGL first, ECRAN second

Data_fem %>%
  select(ID, OXT, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = OXT, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - ECRAN first, OGL second")

## Females - ECRAN first, OGL second

Data_fem %>%
  select(ID, OXT, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = OXT, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

7 Wilcoxon paired on Order of Conditions - Stres

cat("## Whole Sample - OGL first, ECRAN second")

## Whole Sample - OGL first, ECRAN second

Data_merged_long_all %>%
  select(ID, Vas_Stres, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Stres, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Whole Sample - ECRAN first, OGL second")

## Whole Sample - ECRAN first, OGL second

Data_merged_long_all %>%
  select(ID, Vas_Stres, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Stres, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - OGL first, ECRAN second")

## Females - OGL first, ECRAN second

Data_fem %>%
  select(ID, Vas_Stres, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Stres, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - ECRAN first, OGL second")

## Females - ECRAN first, OGL second

Data_fem %>%
  select(ID, Vas_Stres, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Stres, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

8 Wilcoxon paired on Order of Conditions - WB

cat("## Whole Sample - OGL first, ECRAN second")

## Whole Sample - OGL first, ECRAN second

Data_merged_long_all %>%
  select(ID, Vas_Bine, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Bine, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Whole Sample - ECRAN first, OGL second")

## Whole Sample - ECRAN first, OGL second

Data_merged_long_all %>%
  select(ID, Vas_Bine, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Bine, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - OGL first, ECRAN second")

## Females - OGL first, ECRAN second

Data_fem %>%
  select(ID, Vas_Bine, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "OGL" & Nr_zi == "1") | (OglEcran == "ECRAN" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Bine, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

cat("## Females - ECRAN first, OGL second")

## Females - ECRAN first, OGL second

Data_fem %>%
  select(ID, Vas_Bine, OglEcran, PrePost, Nr_zi) %>%
  filter((OglEcran == "ECRAN" & Nr_zi == "1") | (OglEcran == "OGL" & Nr_zi == "2")) %>%
  drop_na() %>%
  group_by(ID) %>% 
  filter(n() > 3) %>%
  ungroup() %>%
  wilcoxon_paired_func(id_var = ID, time_var = PrePost, value_var = Vas_Bine, assum_check = TRUE)

 Summary statistics 

 Wilcoxon signed rank test on paired samples 

 Effect size 

[[1]]

[[2]]

9 Session Info

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

locale:
[1] LC_COLLATE=Romanian_Romania.1250  LC_CTYPE=Romanian_Romania.1250    LC_MONETARY=Romanian_Romania.1250 LC_NUMERIC=C                     
[5] LC_TIME=Romanian_Romania.1250    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rlang_0.4.6        rstatix_0.5.0      GGally_1.4.0       rio_0.5.16         plyr_1.8.6         summarytools_0.8.8 DT_0.5             ggpubr_0.2.5      
 [9] magrittr_1.5       broom_0.5.6        papaja_0.1.0.9842  psych_1.9.12.31    forcats_0.5.0      stringr_1.4.0      dplyr_0.8.5        purrr_0.3.3       
[17] readr_1.3.1        tidyr_1.0.2        tibble_3.0.0       ggplot2_3.3.0      tidyverse_1.3.0    pacman_0.5.1      

loaded via a namespace (and not attached):
 [1] nlme_3.1-140       bitops_1.0-6       matrixStats_0.54.0 fs_1.4.1           lubridate_1.7.4    RColorBrewer_1.1-2 httr_1.4.1         ggsci_2.9         
 [9] tools_3.6.1        backports_1.1.6    R6_2.4.1           DBI_1.0.0          colorspace_1.4-1   withr_2.1.2        tidyselect_1.0.0   mnormt_1.5-6      
[17] curl_4.3           compiler_3.6.1     cli_2.0.2          rvest_0.3.5        xml2_1.3.1         sandwich_2.5-0     labeling_0.3       scales_1.1.0      
[25] mvtnorm_1.1-0      digest_0.6.25      foreign_0.8-71     pkgconfig_2.0.3    htmltools_0.4.0    dbplyr_1.4.3       htmlwidgets_1.5.1  readxl_1.3.1      
[33] rstudioapi_0.11    pryr_0.1.4         shiny_1.2.0        farver_2.0.3       generics_0.0.2     zoo_1.8-4          jsonlite_1.6.1     crosstalk_1.0.0   
[41] zip_1.0.0          car_3.0-7          RCurl_1.95-4.11    modeltools_0.2-22  rapportools_1.0    Matrix_1.2-17      Rcpp_1.0.4.6       munsell_0.5.0     
[49] fansi_0.4.1        abind_1.4-5        lifecycle_0.2.0    multcomp_1.4-8     yaml_2.2.1         stringi_1.4.6      carData_3.0-2      MASS_7.3-51.4     
[57] grid_3.6.1         parallel_3.6.1     promises_1.0.1     crayon_1.3.4       lattice_0.20-38    splines_3.6.1      haven_2.2.0        pander_0.6.3      
[65] hms_0.5.3          knitr_1.28         pillar_1.4.3       ggsignif_0.4.0     stats4_3.6.1       codetools_0.2-16   reprex_0.3.0       glue_1.4.0        
[73] data.table_1.12.8  modelr_0.1.6       vctrs_0.2.4        httpuv_1.4.5       cellranger_1.1.0   gtable_0.3.0       reshape_0.8.8      assertthat_0.2.1  
[81] xfun_0.13          openxlsx_4.1.0     coin_1.2-2         mime_0.9           xtable_1.8-4       later_0.7.5        survival_2.44-1.1  TH.data_1.0-9     
[89] ellipsis_0.3.0    

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com